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ABSTRACT 



A finite element method for structural optimizatic:'. of 
a prismatic beam of a homogeneous, isotropic material is 
developed. The beam has a rectangular cross-section c: 
constant fixed height, a fixed length, and a fixed volume. 
Structural optimum is defined as that shape which allc-s a 
maximum load vjithin the elastic range. 

A com.puter program is developed to solve the resulting 
system of equations and various example problems are s rived. 
Comparison is made with exact optimum beam designs v/here 
possible. 

The finite element model is able to solve problem.s with 
any boundary conditions and types of loading that are :on- 
sistent v/ith the number of elements selected. 
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INTRODUCTION 



I . 

Optimization has a history as old as the universe itself. 
The planets have taken their optimum position in the solar 
system, the solar system in the galaxy, etc. The whole 
history of evolution is one of gradual optimization. Man 
has used optimization in one form or another since the 
earliest of times. Individual man was weak with respect to 
his hostile environment. He learned to become a social 
animal so that as a group he was strong - optimization. 

In wars between groups spears defeated rocks , arrows defeated 
spears, and so on to the present day. All were essentially 
optimization . 

The early Egyptians formalized some optimization tech- 
niques in their development of geometry and trigonometry. 

The development of the Calculus enabled us to find optimum 
function values. The Calculus of Variations made possible 
the optimization of functionals, the optimum function from 
a function space. 

In structural optimization, however, oftentimes the 
Calculus of Variations produces nonlinear differential 
equations that cannot be solved in closed form and the 
numerical methods of solution are difficult. If these 
differential equations could be transform.ed into algebraic 
equations, a solution might be more readily obtained. 
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The finite element method has been applied to the 
theories of elasticity [1] , [2] and mechanical vibrations 

[3]> as well as to other fields, to resolve such problems. 

It has been used in the optimization of truss networks but 
as yet not in actual beam optimization. Since the m.ethod 
is essentially variational in nature, it should have an 
application to structural optimization theory. This sug- 
gests the desirability for research in this field. 

Since this is a first attempt at a formulation of this 
type, we will restrict ourselves to a specific type of 
problem; the structural optimization of a prismatic beam 
of one material which is homogeneous and isotropic. This 
beam has a rectangular cross-section of uniform fixed height, 
and has a given volume and length. The structural optimum 
in this case is the shape required to maximize the load 
within the elastic range. 

A variational formulation is first presented to establish 
the optimization technique for isoperlmetrlc problems. The 
finite element method is applied to the problem using a 
parallel approach. A model is developed and then checked 
for validity by solving various applicable problems. 



10 



II. VARIATIONAL FORMULATION OF BEAM OPTIMIZATION PROBLEMS 



A. DEFINITION OF THE PROBLEM 

Salinas [^] shows that the problem of structural opti- 
mization of a beam, under a fixed volume constraint, may be 
solved by using calculus of variations techniques. He shows 
that operation with the techniques for isoperimetric prob- 
lems [ 5 ] on the augmented potential energy functional, 

T* = T-XS, produces equations which describe the struc- 
turally optimum beam. T is the potential energy of the beam 
and S is the isoperimetric constraint. The quantity X is 
a constant Lagrange multiplier which is actually a measure 
of the strength of the beam. 

We shall consider the problem of structural optimization 
of a simply supported beam of length, L, and volume, . 

The applied load is uniform, P^ (force/unit length). The 
cross-section of the beam is such that the moment of 
inertia may be given as I = Ca'^ where C and n are constants. 
Gravity forces and shear effects are assumed negligible and 
the usual beam assumptions [6] are used, 
i) small deflection theory 
il) uniaxial stress state 

Figure 1 graphically depicts the beam and defines the x 
and y axes. We seek the shape of the beam, to give a maximum 
load intensity subject to the constant volume constraint. 
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Figure 1. Simply Supported Beam, 



B. THE VARIATIONAL FORMULATION 

We first form the potential energy functional from 
strength of materials considerations. The displacement of 
the beam Is denoted by v and the primes Indicate differen- 
tiation with respect to x. 



T 



L 

/ 

0 




P V 
o 



dx 



L 

The fixed volume constraint m.ay be written as / A(x)dx=V 

0 ° 

thus the constraint augmentation becomes: 



L 

As = A/ A(x) dx 
0 



Noting that I(x) may be denoted as CA(x)^, the augmented 
potential energy functional Is: 



T« = / 

0 L 



ECA(x) 
2 



n 



(v")^ - P^v - Aa(x) 



dx 



where v and v" are v(x) and v"(x) the displacement and 
curvature of the beam respectively. 
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In order to extremize the above functional, v;e first 
take the variation [7] of the functional v;lth respect to 
displacement, v, and set this equal to zero. The result 
is the differential equation of equilibrium. This variation 
yields : 

(EAC^v")" - = 0 

which may be written: 

(A^v")" = Pq/EC . 



Integrating this equation twice with respect to x gives : 

/ 2 \ 

A^v" = P /EX — + K^x + K. 

o \ 2 1 2 

The constants of integration, and K^, are obtained 
from the boundary conditions. For a simply supported beam, 
the end moments are zero, or: 

(a) A^(0)v"(0) = 0 

(b) a’^(L)v"(L) = 0 

Boundary conditions (a) and (b) yield K 2 = 0, = -L/2 

respectively . 

Thus we obtain the differential equation of equilibrium 
for the simply supported beam. 



A^v" 



P 

o 

2EC 



( x^-Lx ) 



(II-l) 



Taking the variation of the functional with respect to 
area, A(x), and setting this equal to zero, gives the dif- 
ferential equation which we refer to as the optimality 
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condition, For our exarple the optirallty condition is: 



^ (A)^"^(v")^ - A = 0 

Solving for v" yields: 

_(nzi) 

v" = / 2A A ^ ' (II-2) 

VnEC 

n-1 

Equation (II-2a) shows that, for an optimum beam, the 
curvature is constant only for n = 1. Otherwise the product 
of a power of the cross-sectional area times the curvature 
is constant. 

Multiplying both sides of equation (II-2) by A^ and sub- 
stituting into equation (II-l), the development continues as: 



a'^v" = 



From equation 



n+1 




(II-l): 



n+1 

A 2 flK 
" v/nEX 



P 

o 

2EC 



(x^-Lx) 



Hence , 



(II-3) 



,n+l 



nP 



Be^ 



° (x^-Lx)^ 



2 2 2 2 

Since (x -Lx) = (Lx-x ) it is convenient in further 
development to use the latter expression. Substituting this 
into the previous equation yellds: 



1^4 



A 



n+1 



n; 



O /T 2s‘ 

(Lx-x ) 



8ECA 

Taking the (n+1)^^ root of both sides: 
A = 



8eca 



o 2.2 

(Lx-x ) 



n+1 



Introducing the constraint equation and solving for A we 
obtain : 

1 



L L 

/ A(x)dx = / 

0 0 



■"nP^ 



8eca 



° (Lx-x"^) 



n+1 



dx = V 



o 



But since A is a constant we obtain, 



A = 



1 










1 


1 


,Bec 


V 




1 / 


o i 



2 



n+1 



p n+1 

/ (Lx-x*^) dx 



(II-5) 



Rearranging (II-5) and solving for we obtain the 
equation defining the maximum load Intensity 

n+1 



P = 
o 



J- 



8ECA , . 2 

n o 



/* f T P.n+l j 
/ (Lx-x ) dx 



n+1 



(II-6) 



Substituting (II-5) Into we obtain for the opti- 

mum area: 

2 



A = 



■TT / T 2 \ n+ 1 

V (Lx-x ) 
o 



(II-7) 



/ (Lx-xL"’'Lx 
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It should be noted here that the integral is not defined 



in closed form for all n and therefore must be numerically 
approximated for som.e problems. 

Summarizing the equations that define the solution to a 
given problem we have: 

V 

A = — 2 ^ (II-8 

L „ 

/ (Lx-x^)^ "^dx 
0 



P 

o 



n+1 

/HcT ( , 2 

v n o 



/ (Lx-x^ )'^^^dx 
0 



n+1 

2 



(II-9) 



The value of X is determined as follows. Fromi equation 
(II-2a) : 



v" = 



2X 



nECA 



n-1 



Prom strength of m.aterials theory: 



M 

max 



E(Iv") 



max 



If we define the yield stress of the beam material, S , as 
the maximum allowable stress, we obtain: 

1m| R 

c, ' max m.ax 



where R is the vertical distance to the furthest fiber 

IT13, X 

from the neutral axis . Thus : 



l6 



E(Iv") R 

S = rax_r^ ^ ^ 

y I max r.ax 



Substituting for v" : 



S = ER 



2X 



m.ax 



EGA 



n-1 

‘max 



S = 



2AE 

n 



R 



/max 



,2 

2E 






R 



(II-IO) 



y 



max 



CA^~1 

It may be shown by example that the quantity ^ Is 

R^ 

Independent of x and thus that A Is Indeed a constant for a 
given type of cross-section. 

The first example Is a constant height rectangular cross- 
section of variable width. The height and width are h, and 
b(x) respectively. 



I(x) = = b(x)h • ^ 

A(x) = b(x)h 

.2 

I(x) = ^ A(x) 



C=^, n=l, R=h/2 

CA^~^ = _C_ = h^/12 ^ 1 
R^ r 2 h^/4 ^ 
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Thus for a rectangular cross-section as above, equation 



( II-IO ) becomes : 



X 




(II-ll) 



The second example is a circular cross-section of 
radius r(x) . 

Kx) = = (.r(x)2)" ■ (i) 

2 

A(x) = 1Tr(x) 

CA^~^ = CA ^ IJ^\ 

p2 I^tt/ 



= r 



(wr^) 



1 

ij 



Thus for a circular cross-section, equation (II-IO 
becomes : 



X 




(11-12) 



The third example is a rectangular cross-section of con- 
stant width, b, and variable height, h(x). 

Kx) = =: t,3h(x)3 • ^ 

12b 

A(x) = b • h(x) 

I(x) = — ^ A(x)^ 

12b 
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c = 



12b ‘ 



n = 3, P = h/2 



CA 



n-1 



R 



CA‘ 

r2 



12b' 



(bV) 



1 

3 



Thus for a rectangular cross-section of constant width 
and variable height, equation (II-IO) becomes: 

o 2 

(11-13) 



X = 



-2L 

2E 



Note that equations (II-ll), (11-12), (11-13) depend 

only on cross section type and not on the manner of loading 

L 

(Independent of the term / P(x)v(x)dx). They are also 

0 

Independent of the boundary conditions since they were 
derived from equations in which the boundary conditions had 
not been envoked. This observation v■^ill be used in Section 
III. 

Equations (II-8), (II-9), and (TI-10), when applied to 
a specific problem, yield the structurally optimum beam 
shape for maximum load Intensity for beams which conform to 
the assumptions made at the beginning of this section. These 
beams must be simply supported and under a uniform load 
distribution . 

Appendix A gives three examples of structural optimiza- 
tion problems solved with equations (II-8, 9, 10). 
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III. A FINITE ELE!'-''ENT MODEL FOR BEAM OPTIMIZATION 



A. SPECIFICATIONS AND ASSUMPTIONS 

This development parallels that of the continuous vari- 
ational formulation of the preceding section. V^e consider 
prismatic beams of one material which is homogeneous and 
isotropic. The modulus of elasticity and yield stress of 
the material are E and S respectively. The cross-sectional 
area of the beam is rectangular with constant height, h, and 
variable width, b(x). The length and volume of the beam are 
denoted by L and respectively. 

The plane of loading is in the centroldal plane parallel 
to the y axis. Gravity forces- and shear effects are assumed 
to be negligible. The usual beam assumptions [6] are taken: 
i) small deflection theory 
il) uniaxial stress state 

B. STATEMENT OF THE PROPLEM 

As in Section II, we seek to optimize the shape of the 
beam (cross-sectional area as a function of position along 
the beam axis) to give a maximum load intensity subject to 
the constraint that the beam volume is fixed. The finite 
element method results in an approximate solution to this 
problem. 
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C. DEFIKITIOn OF Ti:~ FINITE ELEFF 
The bean, of Fip^ure 1, vjithout supports, is discretized 

into finite elements of equal length as shown in Figure 2. 

The elements are Identified by numbering them sequentially 
from the origin. The x (centroldal) axis and y axis (posi- 
tive dov;nv;ard) define the plane of loading and bending. 

V/e consider a typical element, say the i^^, as shovjn in 
Figure 3* The displacements and slopes at the left and 
right ends of the element are denoted by 

respectively. The cross-sectional area is . The local 
coordinate system, x. and y. are defined for each element. 

D. APPLICATION OF THE FINITE ELEMENT METHOD 

In accordance with the finite element method [1] , shape 
functions are assumed for the displacement of each element. 
Since a cubic polynomial function for v^ , the displacement 
over the i^^ element, yields continuity of displacements 
and slopes at each elem.ent boundary, it is taken as the 
assumed displacement field. Thus the displacement over the 
i^^ element may be written; 

v(x. ) = <N(x. )> {v} . 

1 1 1 

1x4 4x1 

where <N(x^)> is the cubic shape function vector as developed 
in Appendix B-1 and {v}^ is the column vector (v^, 6^, ^2,02^ 
The elemental cross-sectional area. A., is assumed to be 
constant over the element to obtain a tractable mathematical 
model. This point is further discussed later in the 
development. 
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Figure 2. A Finite Element Beam. 




Figure 3 . 



T-th 

The I 



Finite Flement. 
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Since the augmented potential energy functional, 

= U-V/-AS, must be extrem.ized, finite element formulations 
must be obtained for the: 

i) strain energy of bending (U) 

ii) potential energy of external force (W) 

iii) Isoperimetrlc constraint ( S) 

1 . Strain Energy of Bending 

The strain energy of bending of a beam may be given 

as : 

U = / ^ (v")^ dx 

0 



The finite element assumption that each element is 
under uniaxial stress enables the strain energy in the i^^ 
element, u^ , to be expressed as: 



u . 

1 



A El. 

: V 




dx. 

1 



where vV denotes the second derivative of v(x. ) with respect 

to X. . 

1 

For a rectangular section of constant height, h, and 
width, b^, the moment of inertia of any section of the 
element may be written as : 

= CA^ C = h^/12 

The strain energy in the i^^ element is obtained as 

follows : 

v(x. ) = <N(x. )> {v} . 

1 1 1 
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v"(x^) = <!j"(x^)> 

(v")^ = {v}T <M">T <M">. {v}. 

1 1 11 



u 



1 



/ ^ A 

0 2 *1 



{v}T <N">T <M">. {v}. dx. 
1 1 111 

1x4 4x1 1x4 4x1 



Since the product of the second derivative of the 

T 

shape function vector, <NV NV>, is the only function of 

x^ in the expression, the remaining terms may be brought 
outside the integral, and the integral may be evaluated. 
Thus . 



u . 
1 




{v}: 



/ 

0 



<N">: <N">. dx. {v} 



The Integral is evaluated in Appendix P-2 and called [k^] , 
a 4x4 matrix. Therefore: 



“i ” T *1 ‘''’i '''h 

1x4 4x4 4x1 

Here it may be noted that, if the has a shape function 
vector associated v/lth it other than a constant, the vector 
{v}^ could not be removed from the integrand. Since the 
values for the vector are unknown, this vector would 

appear in the Integrand between two shape function vectors, 
and therefore closed form integration is not possible. 

When the term ECA^ is Included in the Integrand, the 
result of integration is commonly known as the elem.ent 

In the present formulation [k^] is called 
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stiffness matrix. 










iiiMtai 











are 



the "rrodified elerent stiffness r.atrix" since the kv . 

^ J 

measures of elem.ent stiffness. 

Since strain energy Is a scalar, the total strain 
energy In the beam is the sum of the Individual element 
strain energies. Hence we have: 



U = 



n 

E 

1 = 1 



u. = 



n 

E 

1 = 1 



^ A 
2 *1 



{v},: [k*] {v}. 



where n is the number of elements. 

2 . Potential Energy of External Forces 

The potential energy of external forces for the 1^^ 
element may be written as : 



w , 



P <D>. {v}. 
oil 

lxi| 4x1 



The vector <D.> Is called the element consistent load dls- 
1 

tribution vector and Is described In detail In Appendix B-3- 
Since the potential energy Is also a scalar, the 
total Is the sum of the contributions of the individual 
elements. Therefore: 



W = P 



n 

E 

1 = 1 



<D>. {v}. 



3 . The Augmented Potential Energy Functional 

The constant volum.e constraint In finite element 
n 



form Is 
sity P 

o 

n 

T« = E 
1 = 1 



E A.A = V . In order to maximize the load Inten- 

1 = 1 1 o 

we form the augmented potential energy functional. 



Y ''i <''>1 



(V), - 



n 

E 

1 = 1 



o 



<D>. {v}. 
1 1 



n 

X E 
1 = 1 



A.A 

1 
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E. THE FINITE ELEMENT VARIATIONAL PROPLEF 



The extremization of the ootentnal energy functional 
must be done globally, l.e., over the entire beam, A trans- 
formation is therefore made from the local coordinates 
to a global coordinate system, {v^} . The transformation is 
made by numbering the components of sequentially from 

the origin, v;here the first component is the dlsplacem.ent 
at the origin, the second as the slope at the origin, etc. 
Thus ; 




Therefore the global representation of the local "dis- 
placement" vector becomes: 

'''h ' p21-l] = 

^2i + l 
^^ 21+2 



In global coordinates the potential energy functional 
becomes ; 



T» = 



2 



n 

Z 

1 = 1 



A. 

1 



{v.}^ 

1 



[k«] 



<vp - 



n 

Z 

1=1 



<D> 



up 



n 

AA Z A. 
i = l ^ 
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The algebraic equations describing th-- structurally 
optimum beam shape and maximum permissible load intensity 
arise from taking the partial derivatives of the global 
augmented potential energy functional with respect to the 
state variables of individual slopes and displacem.ents and 
also with respect to the control variables of the elem.ent 
areas and setting these equal to zero. The constraint 
equation is also imposed. 

The following set of equations, as developed in 
Appendix B-4, results. The equations must be numbered as 
shown to facilitate enforcement of boundary conditions . 
This is described later in the development. 




( 1 ) 




( 2 ) 






( 1 - 1 ) 



EC 




2jMj+i-2) 



(i) 




^o^(2n+l) ^ 
EC 



0 



(2n+l) 
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0 



(2n+2) 



T — 



•n jf, ^(j+2n-2) 



^o^(2n+2) 

EC 



4 4 

2 2 kf . V / . 



1 = 1 . = i"ij ''(j + 2m-2)^(l+2m-2) 



2XA 



= 0 



m = 



EC 

1,2,. . . ,n 



(2n+2+m) 



n V 

2 A - = 0 (3n+3) 

1 = 1 1 A 



where : 



D. 

1 

A 



V . 

1 



1_ ^ 

1 component of the 
distribution vector, 



global consistent 
<n > 

lx(2n+2) 



load 



strength constant S /6E as developed in 
Section II. ^ 



t h 

i^ component of the global displacement 
vector . 



This is a set of 3n+3 nonlinear algebraic equations 
which, when solved simultaneously, results in: 

1) values for n+1 nodal displacements and n+1 nodal 
slopes at maximum load Intensity. 

11) the n optimum element areas 
ill) the maximum load intensity, 

No boundary conditions are included in the preceding 
global formulation. Since the derivatives may only be taken 
with respect to variables , a fixed displacement or slope 
boundary condition requires the removal of the equation 
associated vjlth the variation of T* with respect to that 
boundary condition. Its value must be substituted In the 
remaining equations. Boundary conditions must be consistent 
with the finite element model, l.e. they must occur at 
nodal points. 
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For example. If the state variable, v. , is fixed, 

the equation is eliminated from, the set, and the value 

for V. substituted in the remaining equations. 

J 

In this way, boundary conditions are Included in the 
formulation for specific problems. 

The set of equations with boundary conditions considered 
is the finite element model for the approximate structural 
optimization of a particular rectangular cross-section beam 
of uniform height and fixed volume. 
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IV. APPLICATION OF THE FINITE ELEMENT MODEL 



In order to apply the model developed in the preceding 
section we must find some method for solving the system of 
nonlinear simultaneous algebraic equations. The size of the 
system for even a small number of elements requires the use 
of a computer algorithm for solution. The only such 
algorithm available to the author was a Fortran IV-G sub- 
routine called SUBROUTINE NLNSYS. This subroutine is 
included in the computer program section of this report. 

A. SUBROUTINE NLNSYS 

This subroutine, in theory, solves simultaneously, by 
an iterative scheme developed by Brown [8], [$] , [10], a 

system of nonlinear algebraic equations of any order. As 
it is dimensioned it is restricted to, at most, 100 
equations . 

NLNSYS requires as input a vector of starting values 
for all unknowns. This starting guess vector must be in 
the region necessary for convergence or the process will 
diverge. This region cannot be defined and different start- 
ing guesses may be necessary to obtain convergence. It 
should be noted that large systems of equations require 
large starting guess vectors. Since the probability of 
a component being outside the region of convergence Increases 
with the number of components, the probability of a starting 
guess vector being within this range decreases with 
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increasing system size. This v:as found to be a difficulty 
that limited the size of the example problems to three 
elements . 

NLNSYS also requires an external subroutine that speci- 
fied the system of equations. Because the system of 
equations developed in Section III changes with the number 
of elements and boundary conditions of a problem, a general 
subroutine is desired that v;lll generate the proper set of 
equations for a given problem. SUBROUTINE FCMLST , provided 
in the program section, accomplishes this task. 

B. SUBROUTINE FCNLST 

SUBROUTINE FCNLST generate a set of equations which are 
evaluated at a point provided by NLNSYS. NLNSYS selects 
only one equation at a time, thus the set of equations must 
be ordered in some manner. 

SUBROUTINE FCNLST sequences the equations for a given 
number of elements as they are ordered in Section III. It 
evaluates these equations at the point specified by NLNSYS. 

The equations derived from parameters that are fixed 
by boundary conditions are removed from the general set and 
NLNSYS operates on one of the equations in this reduced set. 

Problem parameters are read into the program in SUBROUTINE 
FCNLST the first time that it is called by NLNSYS. 

C . SUBROUTINE OPBEAM 

SUBROUTINE OPBEAM, as provided in the program section, 
reads the initial guess vector into the program and provides 
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the output from the program. It permits an efficient m.eans 
of trying multiple starting guess vectors for a given 
problem. It starts the solution to the problem by calling 
NLNSYS. 

D. INSTRUCTIONS FOR PROGRAM USE 

It is first necessary to define the problem to be 
solved. ¥e specify: 

1) the length of the beam (in) - BLENG 

2) the volume (In^) - VOL 

3) the height (in) - HI 

k) the modulus of elasticity (Ib/ln^) - E 

5) the yield stress (Ib/in^) - SIGYP 
Selecting the number of elements, NI , we number the dis- 
placements and slopes at element boundaries. This must be 
done as Indicated in Figure 




Figure H. The Finite Element Numbering Scheme. 



32 



We consider the boundary conditions of the problein and 
form two vectors. One, called IPC, is a vector of subscripts 
of the components of the displacement/slope vector, V, which 
are fixed. These are ordered in Increasing value. The 
second, called BC , is a vector of the values of these fixed 
components . The number of boundary conditions is called 
NBC. For example, if the displacements at the left and 
right ends of the beam are fixed at zero and the slopes at 
these points are variable (simply supported), the vectors 
are : 

IBC(l) = 1 BC(1) = 0 

IBC(2) = (2 X NI) + 1 BC(2) = 0 

The total number of parameters is equal to (3 x NI) + 3 
and is denoted NT. Subtracting the number of boundary con- 
ditions, NBC, we obtain the total number of unknowns, 

NU = NT - NBC. This is also the size of the reduced system. 

Prom the given loading condition on the beam we form the 
global consistent load distribution vector <D> for the number 
of elements chosen. This is denoted CLOAD in the program. 

The initial guess vector, X, is formed by guessing the 
magnitude of the unknown displacements and slopes, the 
areas of the elements, and the maximum load intensity. 

These must be ordered as follov;s : 

1) Displacements and slopes as ordered in the 
vector {v} v;lth boundary conditions removed. 

2) Element areas starting at the origin. 

3) Maximum load intensity. 
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We select values for WAXIT, NUWflG, and IPRIIIT fror. 
their descriptions in SUER0UTI!1E NLffvc. 

A calling program is v;ritten which dimensions X 
(the starting guess vector). This dimension is NU . The 
program calls SUBROUTINE OPPEAM with num.erlcal values for 
the arguments except for X. A sample calling program is: 
DIMENSION X(NU) 

* 5 ^ ^ ^ 

CALL OPBEAM ( NU ,MAXIT , NUMSIG , IPRINT ,X ) 

STOP 

END 

* (numerical values) 

The user must provide the two statements for SUBROUTINE 
FCNLST as described in the program. The first is the 
dimension statement: 

DIMENSION IBC(NBC) ,BC(NBC) ,CLOAD( 2xNI+2 ) ,G (NT) , 
H(NT) ,V(2xNI+2) ,A(NI) 

The dimensions are numerical values. 

The second is the data statement: 

DATA NU/«/,NBC/*/,NI/*/,NT/*/ 
where * indicates the numerical values for the parameters. 

Data cards are prepared and inserted in the following 
order and format : 

1) The guess vector X - 6E12.5 

2) The physical constants: 

BLENG,E,HI,VOL,SIGYP - 5E12.5 
The boundary conditions: 
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IBC(l) , BC(1) - T10,E12.5 
IBC(2) , BC(2) - I10,E12. 5 



IBC(NBC) ,BC(NBC) - I10,E12.5 
4) The consistent load distribution vector 
CLOAD - 6E12.5 

The final output may be one of the following: 

1) A singularity was generated at the output vector. 

2) The Iteration limit was reached. 

3) A solution v;as reached. 



In case 1 a new starting guess should be made and the 
problem re-run. In case 2, if the number of iterations was 
small (<30)j try a larger number of Iterations. If it was 
not small, try a new starting guess. Many guesses may be 
required before a solution is reached. 

An example problem set-up is shown in Appendix D along 
with the results obtained for different problems. 
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V. RESULTS, CONCLUSIONS, AMD RECOI-'I^^ENDATIO^ts 



A . RESULTS 

Figures 22 , 26 , and 28 compare the three element finite 
element optimum shape with the exact optlm.um shape for a 
specific beam under a specific loading condition. Table I 
compares the maximum load intensities for the finite element 
designs of Figures 22-28 with those for the beam of equal 
volume, length, and height but constant cross-section. 



TABLE I 

COMPARISON OF MAXIMUM LOAD RESULTS 



Loading 


A 




B 


c 


D 


E 


Simply Supported 
Uniform Load 


510 Ib/ln 


764 


Ib/ln 


721 


Ib/ln 


l4l 


5.6 


Concentrated 
Loads at 
















X = 10 in 


100 KIP 




— 


194 


KIP 


194 


— 


X = 20 in 


55 KIP 




- 


105 


KIP 


191 


- 


X = 30 in 


40.7 KIP 




- 


76 


KIP 


l'87 


- 


X = 40 in 


34.4 KIP 




- 


63 


KIP 


183 


- 


X = 50 in 


31. 4 KIP 






58 


KIP 


181 




X = 60 in 


30.6 KIP 




- 


57 


KIP 


186 


6.6 


Triangular Load 


996 Ib/ln 




- 


1830 


KIP 


185 


- 


Cantilever Beam 
Uniform Load 


191 Ib/ln 


382 


Ib/in 


360 


Ib/ln 


188 


5.6 



A = Maximum Load for Constant Area Beam 

B = Maximum Load for Exact Optimum Beam 

C = Maximum Load for Three Element Beam 

D = % Increase of Load from Case A to Case C 

E ■= % Deviation of Case C from Case B 
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It is clearly evident in Figures 22 , 2c, an 1* that the 
three-elenent design approximates the exact opti: shape. 

The best criterion for the accuracy of the apprc xi~.ation , 
however, is not the beam's appearance. Since our roal is to 
maximize the load intensity, the accuracy of the approxlm.a- 
tion is best measured by the closeness of the Iced intensity 
allowed by the finite element design to that alli'red by the 
exact optlm.um design. 

In those examples where exact solutions were known, even 
the coarse three-element approximation allows a maximum load 
intensity that is close to that of the exact optimum shape. 
Table I shows that the three-element design gives maximum 
loads within 5.6% of the exact design maximum fcr the simply 
supported and cantilever beams under uniform load and within 
6.6% for the simply supported beam under a concentrated 
load at its center. 

Attempts were made to Increase the number cf elements 
in the simply supported beam under a uniform lead. This 
was desired in order to show that the shape would converge 
to the exact solution with an increasing number of elements. 
However, difficulty was encountered in finding an initial 
guess vector within the range of convergence. lue to the 
limited time available and the need to check the model's 
validity for other loading and boundary conditions, the 
three-element model v/as used throughout. 
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B. COMCLUSIOnS 



The above results show that the finite eler-r.t model 
developed in Section III, for the structural c.i t imlzatlon 
of a homogeneous beam of uniform fixed height, and fixed 
volume. Is valid. Furtherm.ore , the results sh that a 
relatively small number of elements may be use! to design 
a beam that will closely approximate the strength of an 
exact maximum strength design. 

As a practical matter, the construction of a beam with 
the shape of the finite element solution Is mu:-h easier 
than constructing the exact optimum shape. In the exact 
optimum design, points where the cross-sectional area is 
zero must have some material added to account for the shear 
effects that were neglected In the formulation. Much less. 
If any, material must be added to the finite element design 
at the same points. The finite element design would have 
stress concentrations at the element boundaries which are 
not present in the exact optimum design. 

The finite element model provides a means of approximat- 
ing solutions to optimization problems which are Impossible 
to solve In a closed form by classical technloues and 
extremely difficult to solve numerically. This Is due 
mainly to the algebraic nature of the finite element method 
and the differential equation nature of the classical 
methods . 

The model developed was, of necessity, fcr the simplest 
type of beam and loading conditions. The appz'-cach used. 
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however, should be valid, v.'lth some modifications, for the 
more complex cross-sections. It has been shown that the 
finite element techniques can be successfully applied to 
the problems of beam optimization. 

C. RECOMMENDATIONS 

Many more problems for which the developed model applies 
should be solved with various combinations of loading and 
boundary conditions for a parametric study of the beam 
optimization problem. Problems with more elements should be 
solved to show convergence to known exact solutions. 

Although SUBROUTINE NLNSYS provided solutions, the 
algorithm is very sensitive to the starting guess in rela- 
tion to the region of convergence. More than one attempt is 
usually required to obtain a solution to a given problem. 

An algorithm that does not exhibit this sensitivity is 
required to enable m.ore efficient problem solving using 
finite element techniques. At the present time very few 
algorithms for solving systems of nonlinear equations exist. 

Another technique that is possible would be to use some - 
type of direct search or steepest descent algorithm to find 
an optimum solution to the finite element potential energy 
functional Itself. VJlth this type of formulation it may 
be possible to allow variable element areas. A solution 
might be more readily attained as compared with a solution 
to the system of nonlinear equations. Some work was started 
in this direction late in the investigation. The initial 
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results looked promising but time restrictions dictated 
that this approach be terminated in favor of com.pleting 
the presented formulation. 

Models for other cross section types could be developed 
using a similar approach to that presented in Section III. 
Also models with other than equal size elements or for 
beams of more than one material would be of interest. 

In the laboratory, an experiment could be run comparing 
the behavior of an exact optimum beam with that of its 
finite element approximations to see if the results compare 
with those predicted by the models. 

The entire area of the application of finite element 
methods to optimization has hardly been touched. The 
finite element method is a powerful tool and more research 
should be done in its application in this field. 
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APPENDIX A 



EXAMPLES OF BEAM STRUCTURAL OPTIMIZATION PROBLEMS 
USING THE VARIATIONAL FORMULATION 



1. Consider a simply supported beam of rectangular cross 
section with constant height, h, and variable width, b(x) 
The beam is under a uniform load distribution and has 
given volume and length, and L respectively. Find the 
optimum shape, b(x) such that the load Intensity P^ is a 
maximum. The beam has a modulus of elasticity E and yiel 
strength S^ . 

Solution : 

The moment of inertia of the cross-section is: 



I 




Thus C = h^/12 and n 
obtain : 



(bh) 

= 1 . 




• A 



Using equation (II-8) vjlth n = 



V (Lx-x^) 
o 



/ (Lx-x )dx 
0 



Invoking equation (II-IO) to obtain (II-ll) the value of 
becomes : 




Substituting the values for C, X, and n into equation 
(II-9) yields : 



^1 



h? 

V o 



3/ (Lx-x'^ ) dx 

0 

The Integral can be evaluated in closed form. 

/ (Lx-x*^) dx = L^/6 

0 

Thus : 

'‘'"’opt - ) 

v;hich gives: 

o 

The beam shape Is graphically -depicted In Figure 5. The 
maximum load Intensity is: 

2hS V 

P = ^ 

o T 3 

max L 

This shows that, for the given problem, maximum load 
Intensity varies linearly with the height, yield stress, 
and volume, and Inversely with the cube of length. 

2. Consider a simply supported beam of rectangular cross- 
section with a constant width, b, and variable height, h(x) 
The beam Is under a uniform load intensity P and has a 
given volume, V^, and length, L. Find the optimum shape, 
h(x) such that the load Intensity, P^, is a maximum. The 
beam has a modulus of elasticity, E, and yield strength, S 

O' 






Solution : 



The moment of inertia of the cross-section is: 



I = 



bh- 



12 2 
12b 



^ a3 



Thus C = r- and n = 3. Using equation (II-8) with n 



12b‘ 



A = 



V^(Lx-x^)'® 

L _ V 

/ (Lx-x'^ ) ^dx 
0 



Equation (11-13) gives : 



X = S^/2E 

y 



With A, C, and n above, equation (II-9) becomes 



S V ‘ 

P = 



o 3b 



” L 2 ^ ' 

/ (Lx-x^) dx 

- 0 



The Integral may be evaluated and is 



L ^ u _T 2 

r f T ^ \ ^ ^ TtL 

/ (Lx-x ) dx = 



0 



8 



Thus : 



opt 



V 

8V^(Lx-x^)^ 2.55 V (Lx-x^)^ 



Solving for h(x) : 



2.55 V^(Lx-x^) ' 
bL^ 



h(x) ^ 



h(x) Is gra|±ilcally depicted in Figure 6. The maxinun per- 
missible load intensity is: 



Po 



max 



6^S 



. 2 ,. 4 

s bL 



2. ITS V 

y o 

bL^ 



Here Po varies linearly with yield stress, para- 

ina.x 

bolically with volume and Inversely with width and the 
fourth pov;er of length. 



3. Consider a simply supported beam of circular cross- 
section under a uniform load distribution. It has a given 
volume and length, and L respectively. Find the optimum 
radius, r(x) such that the load intensity is a maximum. 

The beam has a modulus of elasticity, E, and yield strength, 
S . 

y 

Solution : 

The moment of inertia of a circular cross-section is: 



I 



4 

irr 

4 



1 



(irr^ 



2 



CA 



n 



Thus C 



4it 



A = 



and n = 2 . 
V^(Lx-x^) 



VJlth 

2/3 



L p 2/3 
/ (Lx-x”^) dx 



n = 



•0 

Equation (11-12) gives the value 
section as : 



2, equation (II-8) becomes 



of X for a circular cross- 
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V/ith A, C, and n as above, equation (3-9) tecor.es: 



P 

o 



S V 
y o 

2/tt 



1 

Fl p 2/3 I 
/ (Lx-x ) dx 

Lo J 



3/2 



The Integral may be evaluated in closed form. However it 
Involves the Gamma Function which must be approximated. 
The integral is, then: 



L p 2/3 
/ (Lx-x ) dx 
0 



.0979 • (L)'^/3 



Thus : 



A 



V (Lx-x'^) 
o 



2/3 



.0979l'^'^^ 




(Lx-x ) 
.308 L^^^ 



2/3 



The beam shape is depicted In Figure 7. The maximum 
load intensity is given by : 



Po 



S V 

..y o 



3/2 



.1067 L^'^^ 



Po is linear in yield strength, proportional to the 
m&x 

7/2 power of volume and inversely proportional to the 7/2 
poiver of length. 

Reference [11] gives optimum beam shapes for example 
problems 1 and 2 derived in a different manner. The shapes 



^15 



shown in Figures 
The ir.aximum load 



5 and 6 are identical to those shapes, 
intensities also correspond. 



^6 
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Flfrure P. OotiF’ur ^ear Fhaoe for a Pectanpular Cross- 
section Pear of Uniforr Feipht unPer a 
Uni^orr Load. 
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hbL 




X 



uniforiTi load distribution 




t 

y 

Figure 6. Optimum Peain Fhaoe for a Pectangular Pross- 
section Peam of Uniforr V'idth under 
Uniforr Load. 
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F’ipure 7- Ontjrur’ Pe^r Pbpoe for a Circular Cross- 
section Pear under a I'niforr Load. 



APPEIIDIX P 



CALCULATIONS FOP THE FINITE ELEMENT METHOD 
FOR BEAM STRUCTURAL OPTIMIZATION 



1 . The Cubic Displacement Function 

V/e wish to find a cubic displacement function vector 

<N(x.)> such that 
1 

v(x.) = <N(x.)> 

1 1 

where v(x^) is a cubic polynomial. 





S 2 

VJe let v(x^) = ax^ + 6x^ + + 6 over the Interval 

0 £ x^ £ A . Boundary conditions are: 
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(1) v(0) = 

(2) v'(0) = 0^ 

(3) v(A) = V 2 

(4) V’ (A) = 02 

Substituting boundary condition (1) in v(x^) gives 6 = . 

Taking the first derivative of v(x^) with respect to 
x^ yields : 

v*(x.) = 3ax? + 28x. + y 

Boundary condition (2) in v'(x^) yields: 

Y = 



Substituting the values for y and 5 into v, and v' and 
invoking boundary conditions (3) and (4) yields: 




6 



3v^ 20^ 3 v2 02 

? r " ^ - T 



In vector notation: 



a 

8 



^a3 .2 




A^ A 



2 _1 
" a2 



>{v)^ 



y=<0100> {v}^ 
6=<1000> {v}^ 



Thus in vector notation: 
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v(x. ) 



3x2 

1 1 



/ x .2 27.2 



+ 1 



\A-^ A 
lx} x?\ 






— + X. 

A 1 



^ 2 
2x; 3x/ 






A‘ 



i_ 

2 



_i 

A 



>{v}^ 



Therefore the shape function vector <N(x^)> is: 



<N(x.)> = 
1 

lxi| 



'2x? 3x? 

1 1 + i'/ 1 



/x? 

1 

a2 



3 a2 

X' 



A‘ 

,2 

"i 



3 2 

X? 2x7 

='1, 



2x} 3x? 



A^ 



The second derivative of this vector v;ith respect to x^ is 



<N(x. ) > = 
1 



12x. 
1 



/6x, 



A‘ 



A 



12x 



i . 6 If^^'i _ 2 

a3 a2 a2 



2 . The Modified Element Stiffness Matrix 

The modified element stiffness matrix is defined as : 

A T 

[k*] = / <N"(x.)>'^ <N"(x.)> dx. 

1 1 1 

4xi| ° i]xl lxi| 



Taking the indicated vector product and integrating yields 



[^•] 



12/A^ 


6/a2 


-12/A^ 


6/A 


6/a2 


Va 


- 6/a2 


2/A 


-12/A^ 


-6/a2 


12/A^ 


-6/A 


6/a2 


2/A 


- 6/a2 


4/A 



3 . The Consistent Load Distribution Vector 

t h 

Assume that the load per unit length over the 1 element 
may be written as an r^'^ order polynomial 
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P(x^) = a + 3x^ + Yx^ + ... + nx^ 

where a, 6, y , . . . ,r\ are constants v/hlch are determined from 

knovm values of load intensity at r+1 loactions In 

0 £ 1. ^ • fhese known values are represented as fractions 

of the maximum, load Intensity, P , over the beam (0£X£L) 

PCXj^) ° 

such as a = a * — where 0 < x. < A. Thus: 

Po - 1 - 

P(x.) = P [a + 3x. + yx? + ... + nx^] 

1 o 1 ' 1 

The load on an increment dx^ of the 1^^ elem.ent is : 

dP(x. ) = P [a + 6x. + yx? + ... + nxT] dx. 

1 o 1 'l 1 1 

The work done on the i^'^ element by this increm.ental force 
is : 

dw . = dP ( X . ) V ( X . ) 

1 1 1 



but : 



v(x. ) = <N(x. )> {v} . 

1 lx 

Ix^ ^^xl 



Thus : 

dw. = P [a+6x.+yx?+ . . .+nxr] <N(x.)> {v}. dx. 

~1 O Dm 3— 3— 3— 3_ 

The work done on the i^^ element is: 



w . 



1 



A 

/ 

0 



dw. 

1 



? r 

w. = / P [a+6x . +y X . + . . . +rix . ] <N(x.)> {v} . dx. 

3_ 3_ 3_ 3_ 3- 3_ 3_ 
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and 



Lettinp [a+ 8x . + yx 7+ . . . +rxf] be denoted bv 0 
noting that and are Independent of , t: 

expression becomes: 

A 

w. = <N(x )> dx {v} 

-L Uq1_L 1 1 1 



The integral is defined as the element consis 

distribution vector <D>. 

1 



A 

<D>. = / Q (x. ) <N(x. )> dx. 
ix^ 0 

. th 



The work done on the i element may then be writ 



w . = P <D> . { v} . 

1 O 11 



The work done on the beam is 



]fi = 



n 

Z 

1=1 



w , 



= P 



n 

Z 

i = l 



<D>. {v}. 



It may be noted there is no requirement for the e 
load distribution polynomial, P(x^), to be the sa 
elem.ent. On the contrary, each element may have 
load distribution polynomial, P(x^), and the disc 
may even be discontinuous at element boundaries, 
requirement is that the ratio of P(local)/P^ be k 
r+1 points over the element in order to define th 
distrlbution polynomial. 

We take the follov/ing case as an example: 

Consider the first three elements of an r. 
beam. Given that the distributions of load on el 



e above 



lent load 



1 en : 



lement 
■"e for each 
a different 
rlbution 

The only 
n own at 
e r^^ order 

element 
em.ents 1, 



2, and 3 are cubic, linear and quadratic respectively, and, 
also given the values of P /P v/here z = A, B, C, 
we seek an expression for the work done on the first three 
elements of the beam. Figure 9 shov/s a graphical represen- 
tation of the problem. 



F 




y ,v 

Figure 9. Multiple Loading Types on a Beam. 

The values P./P through Pr>,/P define the cubic distrl- 
A o Do 

bution function Q 2 (x^) over the first elem.ent. The values 
P„/P and P-c/P define the linear distribution function 

Cj O r O 

Q, (x„) over the second element. The values of P„/P 

1 id r O 

through P^/P^ define the quadratic distribution function 
0,2 (x^) over the third elem.ent. These functions in turn de- 
fine consistent load distribution vectors <D^>, <D 2 > , and 
<D^> for the elem.ents. 

A 

<D,> = / <N(x 2_)> dx^ 

± Q i ± 

A 

<D2> = / Q^(x2) <M(x 2)> dx2 

A 

<D^> = f Q2(y^-^) <N(x 2)> dx^ 
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The v;ork done on the first three elerents is: 



VKl-3) = ^^^3 ^^3^^ 

H . The Concentrated Load 

If the load is concentrated vice distributed, the 
consistent load distribution vector is formed in a manner 
parallel to that of the distributed load. 

We apply a concentrated load of intensity 6 Pq( 0<6£1) 
at a distance aA(0 £«<_!) from the origin of the i^'^ 
element . 




Figure 10. A General Element vjith Concentrated Load. 

The work done on this element is : 

V/. = 3P v(aA) 

10 

But : 

v(aA) = <N(aA)> Wlj, 
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There fore : 



w. = 6P^ <N(aA)> {v}. 

1 o 1 

If we define the element consistent load distribution 
vector as : 

<D>^ = 3 <N(aA)> 

we obtain for the work done on the i^^ element by the load 

P : 

o 



w . 

1 



P <D>. 
o 1 



{v} 



1 



Note that the expression is the same as that for the dis- 
tributed load. However, the units for P and <D> . are 

^ o 1 

different - P is in units of force vice force per unit 
length. 

If multiple concentrated loads are applied to one 
element, the element consistent load vector is the sum of 
those obtained by considering each load separately. In 
symbolic form: 

The individual load vectors are: 

<Di>i = 3^ <N(a^A)> 

<D2>j_ = 3^ <N(k2A)> 

<D.>. = 3. <N(a.A)> 

J J J 

The element consistent load vector is: 

<D>. = Z 3. <N(a.A)> 

1 n* ,1 .-] 
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5 . The Finite Eler:ent Variational Equations. 



\'Je have as the global augr.ented potential energy 
functional : 






n 

P„ E <D>, {v, } 
° 1=1 ^ ^ 



n 

- X E AA. 
1 = 1 ^ 



where n Is the number of elements. 

Define the following as the global consistent load 

distribution vector <D> 

lx(2n+2) 



a 



J . t.2 



<D> = <D^ , D^, + D^, , D^> 



n ^ n ^ n 



— ^ 2 ^ ) D2j*>-jD 



2n+l’ ^2n+2^ 



<D> Is formed by assembling the element distribution vectors 
as Indicated. Thus the global work term becomes: 



n _ 

P E <D.> {v. } = P <D> {v} 

o . , 1 1 o 

lx(2n+2) (2n+2)xl 



The augmented global potential energy functional that will 
be operated on Is: 

= ^ E A. {v.}"^ [k«] {v.} - P^ <D> {v} - XA E A. 

S ^ 1=1 1 1 10 1=1 ^ 

Taking the partial derivative of T* v;lth respect to each 

component v. In {v} and setting result equal to zero gives: 

J 

3T« 

^ = 0; EGA, E k«. v. - P D, = 

1 Ij j o 1 
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0 



0: 



0 



3T« 

9v^ 



EGA. Z k“ . V . - P 

1 j = l 2J J 0 2 



3T» l4 

s— E = 0: ECA £ kS. v,... ,, , 

(1^) 3j (o+i-H) 



+ ECA/./,->\ Z k?. V/.,. — P D/. 

(i/2) Ij (j+i-2) o (i-l) 



= 0 



1 = 4,6,8j...2n 



3T* 4 

5^=0: ECA 



+ EGA/ . Z V, . . . - P D. 

(i/2) 2j (j+i-2) o 1 



1 = 4 ,6 , 8 , . . . 2n 



grp* 



g 


= 0; 


ECA 

n 


Z 

.1=1 


^3j 


^""(2n+l) 


3T* 

g 


= 0; 


ECA 

n 


4 

Z 

J=1 


% 


^""(2n+2) 



o ^(2n+l) 



= 0 



’o ^(2n+2) ° 



This gives 2n+2 equations. 

We next take the partial derivatives of T* with respect 
to the A^'s and set them equal to zero; 



9T 4 4 

-^ = 0: Z Z k¥. V, 



3A. ^ij ''(k+21-2) ""(3+21-2) EC 



2A 



i = 1 , 2 , . . . ,n 



This gives n equations. 

At this stage vie have 2n+2 equations in 3n+3 unknovms 
hence an additional equation is necessary. This equation 



is the volum.e constraint: 



0 



n 

A E 
1 = 

Therefore we have a total of 3n+3 equations in the 3n+3 
unknowns : 

2n+2 components of the vector {v} 

n areas A. 

1 

1 maximum load Intensity 

V/e therefore have a mathematically tractable problem. 
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APPENDIX C 



EXAMPLES OF CONSISTENT LOAD DISTRIPUTION VECTORS 

1. As our first example v;e consider a three element problem 
under a uniform load distribution. 









’ 






i 

1 

! 


2 


3 









o 



-f?* 



^ ■ ' ■■■ '* ■' ( A- «©■ 



A 



Figure 11, A Three Element Beam, Uniform Load. 

Observing the general elem.ent of Figure 12, the work done 
on this element is: 

A 

w. = / P <N(x. )> dx. {v} 



1 Q o 



1 1 




Figure 12. A General Element, Uniform Load. 
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which may be written: 



A 

P / <M(x. )> dx. {v} . 

On 1 1 1 



0 



From our previous definition in Appendix B, the elem.ent 
consistent load distribution vector is: 

A 

<D> . = / <N(x. )> dx. 

Evaluating the integral where <N(x^)> is as defined in 
Appendix B; 



Applying this to the three element problem and noting 
that the load intensity is the same for each element, we 
obtain: 



A 

/ <N(x. )> dx. 

0 ^ ^ 



? 2 
A ^ A _ ^ 

2 ’ 12 ’ 2 ’ 12 



Thus 




(C-1) 





We assemble these element vectors to obtain the global con- 



sistent load distribution vector <D> 



<D> = < 




(C-2) 



The work done on the beam is : 
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/ “ 



Y = 

'' -o 2’ 12 






V, 



V, 



V, 



0^ A3 O3 Q3 



V, 



V/ 



Vr 



V( 



^ ) 



2. As a second example let us take the following problem: 



Po/2 



y .V 





» \ 


• P 

^ 0 


P /2 
0 

r 




^ --1 


1 


2 


3 








»• — 



A — - “A — — A — 



Figure 13 . Three Element Beam, Multiple Uniform Loads 



The element consistent load distribution vector for a uni- 
form load from equation (C-1) Is: 



<D> = < A Al A 

^1 2 ’ 12 ’ 2 ’ 



12 



In this problem, however, the load Intensity Is not the 
same on each element. Therefore the element consistent load 
distribution vectors are: 
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<D>^ = <D> 



= 2 



Or: 



<D>o = <D>. 
2 1 



<D>1 = l> 



2M 



<D> = < ^ ^ ^ 

^2 2 ’ 12 ’ 2 ’ 



12 



Assembling these vectors we obtain for the global consistent 
load distribution vector: 

- A ^ 3A A -A^ 

^ ^ 7T’ T“’ 12’ IP’ 12’ ?’ PP 

3. Let us now look at the following problem with a concen- 
trated load: 




y ,v 



Figure l4. Three Element Beam, One Concentrated Load. 
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lj 










Taking a general eler.ent v;ith a concentrated load: 



Figure 15. 



p 





" o 

» 


1 




• 

A — 




u 

”i ’"i 





General Element , 



X . 

1 



Concentrated Load. 



The consistent load distribution vector for this element is : 



<D>. = <N(aA)> 

1 

But : 

<N(aA)> = < ( 2a^-3ot^+l ) , (a^A-2a^a+aA) , (-2a^+3a^), 
(a^A-a^A) > 

Thus : 

<D>^ = < ( 2a^-3«^+l ) , [ (aA ) (a^-2a+l) ] , [ ( ) ( -2a+ 3 ) ] > 

[(a^A)(a-l)] > (C-il)' 

In the specified problem there are no loads on elements 2 
or 3 thus : 



<D>2 = <D^> = <0, 0, 0, 0> 

The elem.ent consistent load vector for element 1 is there- 
fore : 



<D>^ - <D(-^)>^ 




3A 27 _ 
32 ’ 




(C-4a) 
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Assenbling" obtain the global con- 

sistent load distribution vector <D>. 



<D> 



5_ 3A 22 

32’ 32’ 



QA 

PT’ 



0 , 0 , 



0 , 



0 > 



(C-5) 



A, Let us take example 3 and add an additional 
element 1 of intensity a distance ^ from the 



load on 
left end: 




y ,v 



— ^ X 



Figure l6. Three Element Beam, Tv;o Concentrated 
Loads on Element One . 



The element consistent load 

concentrated load P a distance 

o 



<D( 




< 



27 

32’ 



9A ^ 

P”’ 32’ 



distribution vector for a 
A/^l from the origin is : 




In this example <!)> ^ and <D>^ are 
That is <D >2 = <D>^ = <0,0,0,0>. 
tribution vector for element 1 is 



the same as in example 3- 
The consistent load dls- 
the sum: 



<D>^ = <D(^)> + <D(|^)> 
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Thus : 



<D> 



1 



< 1 , 





> 



Assembling the elem.ent vectors vje obtain for the global . 
vector : 



<D> 





0 , 0 ,0 , 0 > 



(C-7) 



5. For this example let us take the following beam with 
concentrated loads as shown: 




y 



Figure 17. Three Element Beam, Concentrated 
Loads on Elements one and Three. 



From equation (C-6) we have for <D>^ ; 

<D> - < 27 2A 5_ _ ^ > 

^1 32" 6i|" 32" ^ 

Element 2 has no loads thus : 



<D>2 = <0, 0, 0, 0> 



Equation (C-4a) gives the element consistent load vector 
for a loading such as is on element 3 but of intensity . 
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The element consistent load distribution vector for element 
3 Is therefore one-half of that value in (C-4a) . Thus 



<d>2 = < 



3A 27 



9A 



12H’ P"’ IP 



Assembling these to obtain the global consistent load 
distribution vector: 



<D> 



< 



27 9A 5 3A 5 

32 ^ P’ 32" P" P" 



3A 27 9A 

IP" P" IP 



> 



(C-8) 



6 . As a final example let us take a three element beam 
under a triangular load as shown: 



P 




Figure I 8 . Three Element Beam, Triangular Load. 

We first examine a general element with a basic 
triangular load of maximum intensity P as shown: 



68 




!'■ 

■ 1 ' ' 1 



Figure 19. General Element, Triangular Load. 



This distribution may be written: 

Px. 

P(x. ) = 

1 A 

The element consistent load distribution vector by 
definition is: 

A X . 

<D> = f <N(x. )> dx. 

■^0 



Evaluating the Integral we obtain for the basic element 
consistent load distribution for a triangular load, 
P(x^) = Px^: 



<D> . 
1 



3A A^ 21A 

20’ 30’ F0~’ 




(C-9) 



If we depict the example problem as follows: 

p 




y ,v 
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We see that element 1 has a basic triangular load distri- 
bution of intensity Pq/ 3> Element 2 has a uniform load 
distribution of intensity P^/B with the triangular distri- 
bution of element 1 superposed. Element 3 similarly has 
this superposition but the uniform load Intensity is 



2P„/3. 



Prom (C-9) and the relation P = P /3 we obtain: 

o 

<D> = 1 <j5>triang ^ 7A 

<D>i ^ <D>^ - < 20^ 90’ 60’ 20 



From the result for above, and from (C-1) with 



intensity Pq/ 3> the uniform portion of <D >2 is; 
uniform 



<D>: 



2 2 
1 A A -A 

3 2 ’ 12 ’ 2 ’ 12 



Thus ; , 

^j^^uniform 

And : 



2 2 
A ^ -A^ 

S’’ sS”’ F’ 12 



<D> 



2 



<D> 



2 



^P^trlang ^ 
. 13A 7A^ 

I¥o’ 



<D> uniform 

17A 2A^ , 

F0“’ " W~ 



The third element consistent load distribution vector is 
developed in a similar manner; 



And ; 



uniform _ 2^ < A A_£ A ^ 
3 "32’ 12’ 2’ “ 12 

2 2 

_ . A A^ A A^ ^ 

3’ rS"’ 3’ iB" 



<D> = 
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< 



Or : 



<D> 



3 



Assembling <D>^ 
<D> = 



23A a- QA ^ 

^F"’ 15’ 20 ’ 



<D> 2 > and vie obtain: 

^ ^ A A^ 2^ A^ QA 
20’ 90’ 3’ W5’ 3 ’ IT5’ 20’ 
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(C-10) 



APPENDIX D 



FINITE ELEFEMT OPTIMIZATION EXAMPLES 



All of the follov;lng examples consider the approximate 
structural optimization of a rectangular cross-section beam 
of uniform height and fixed volum.e . The physical parameters 
are : 



a) beam length = 120 in 

b) beam height = 6 in 

o 

c) beam volume = 2200 in 

d) modulus of elasticity = 30x10° Ih/ln 

o 2 

e) yield strength = 50x10 Ib/in 

A three element approximation is used throughout. 

1 . Representative problem set-up 

As a representative problem consider a simply sup- 
ported beam under a uniform load. The beam is depicted in 
Figure 21. NI = 3. 




Figure 21. Three Element Problem, Uniform Load. 
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Since the bean Is sinply supported the boundary condi- 



tions are: 



o 

II 

1 — i 
> 




Vy = 0 




Therefore : 




IBC(l) = 1 


BC(1) = 0.0 


IBC(2) = 7 


BC(2) = 0.0 



NBC = 2 

The total number of parameters are: 
NT =3x3+3= 12 
The total number of unknowns are: 

NU = 12 - 2 = 10 



Prom Appendix C, the consistent load distribution vector for 

a uniform load on three elements is: 

2 2 
CLOAD = < ~ > 

Since A = L/3 = ^0 in, CLOAD becom.es: 

< 20 , 133.333, ^0, 0, iJO, 0, 20 , -133-33 > 

The components of the x vector are: 



X = < 



vj^, v^, vg, vg, A^, A^, Ag, > 



We guess values for the x vector: 



= .03 

Vg = .9 in 

Vj^ = .01 

v^ = .9 in 

= -.01 



A^ = 15 In^ 
A^ = 28 in^ 

Ag = 15 in^ 



P = 760 Ib/in 
o 



Vo = 



-.03 
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3 



V.’e select values for 
require 50 iterations 

MAXIT = 50 



"AXI"!" , XUI’flG, 



raxlnun . 



and IPPII'T. 



Let us 



We desire four significant digits. 
NUMSIG = 



We desire the x value printed at each iteration. 
IPRINT = 0 



We now can utilize the program. The calling program is: 

DIMENSION X(I0) 

CALL OPBEAM(I0,50,il,0,X) 

STOP 

END 



We Insert the following cards in their proper locations in 
FCNLST. 

DIMENSION IBC(2) ,BC(2) ,CLOAD(8) ,G(12) ,H(12) , 
V(8),A(3) 

DATA NU/10/, NBC/2/, NI/3/, NT/12/ 

The first data card is the first six components of the 
X vector in E12.5 format. The second is the last four 
components of the x vector in E12.5 form.at . The third card 
is the list of physical parameters BLENG ,E ,HI ,VOL ,SIGYP in 
E12.5 format. The fourth card is the first component of 
IBC and the corresponding com.ponent of BC (I10.,E12 . 5 ) • 

The fifth card is the same as the fourth but for the second 
components of the vectors. 
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m 



The last two cards contain the eig"ht components of the 



consistent 


load distribution vector in 


6 EI 2.5 form.at . 


The 


data deck 


is thus ; 






0. 3E-1 


0.9E0 O.lE-1 0.9E0 


-O.lE-1 


-0. 3E-1 


1 .5E1 


2.8E1 1.5E1 7.6E2 






1. 2E1 


3.0E7 6.0E0 2.2E3 


5.0E4 




1 


0 . OEO 






2 


O.OEO 






2.0E1 


1.33333E2 4.0E1 0. 


OEO 4.0E1 


O.OEO 


2.0E1 


-1.33333E2 







This program was run and the following solution was 
returned after four Iterations . 

= .0310 A = 15.01 in^ P = 720.6 Ib/ln 

^ _L ^ O 

= .9566 in = 24.98 in 

= .0111 = 15.01 In^ 

Vr- = .9566 in 

5 

Vg = -.0111 

Vg = -.0310 

We are Interested in mainly the areas and maxim.um load 
results. If the areas are divided by the height (6 in) we 
obtain the width. Since the beam is symmetric about its 
centerline we may plot one half of the beam to depict its 
shape. Figure 22 is a plot of the three finite element 
approximation v;ith the exact solution superlm.posed . 

2 . Other exam.ple problems 

The problems of a simply supported beam under con- 
centrated loads at various locations were also solved v;lth 



75 




1 



i 




the computer program. The results for these problems are 
shovm in Figures 23, 2H , 25, and 26. Corresponding loads 
to the right of the beam m.idpoint were checked and the beam, 
shapes obtained were the inverse of those shown, i.e., the 
half-width values for the end elements were interchanged. 

The solution of a simply supported beam under a tri- 
angular load distribution is shown in Figure 27. 

Finally the problem of a cantilever beam under a uniform 
load was solved. This is shown in Figure 28 . 

Optimum beam shapes for the simply supported beams under 
1) a uniform load and 2) a concentrated load at the center 
and also for the cantilever beam under a uniform load are 
shown in [11] . These exact optimum shapes are superimposed 
for comparison in Figures 22, -26 , and 28. 
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HALF-WIDTH (in) 




I 1 1 1 — -H— I 1 X 

0.0 20 iin 60 pn inn I2n (in) 




MAXIMUM LOAD, 
MAXIMUM LOAD, 
MAXIMUM LOAD, 



FIHI':^F FLFMENT SOLUTION ... 721 Ib/in 

FXACT SOLUTION Jfl} Ib/ln 

UMIFOPM CROSS-SECTION 510 Ib/ln 



Figure 22. "^hree Flerent Ontirum DesiP’n for a 
Singly supported Pean under 
Uniform Load. 
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FAM HALF-WIDTH (in) 



1 2. 3 4 5 fe 




Loads at points 1 through f' apply in Figures 23 - 26 
The beams shown in these figures have the parameters 
V = 220n in 3 

° A 2 

E = 30x10^ Ib/in 

= 50x10^ Ib/ln^ 



Three Flem'ent Optimum Design for a 
Simply Supported Peam. under, a Con- 
centrated Load Applied at Point 1. 
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Figure 23 . 
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Figure 2^. Three Flenent Optimum Peam. Designs for 
Simply Supported Beams under Concen- 
trated Loads at Points 2 and 3. 
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I 
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Figure 25. Three Eler.ent Optinur Beam Designs for 
Siirply Supported Beans under Concen- 
trated Loads Applied at Points 4 and 5. 
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Figure 26. Three ^lerr’ent Optimum Beam Design for a 
Simiply Supported Beam under a Concen- 
trated Load Applied at Point 6. 
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Flpure 27. Three Element Optimum Hearn Pesipn for a 
Simply Supported Hearn under a Triangular 
Load Distribution. 
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Figure 28. Three Element Optimum Beam Design for 
Cantilever Beam, under a Uniform. Load 
Distribution . 
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